Behavioral Plasticity Analysis

Data preparation

This document includes code to reproduce all analyses and figures in the paper Nielsen, Dablander, Debnath, Emegor, Ghai, Gwozdz, Hahnel, Hofmann, & Bauer (2024). Perceived plasticity of climate-relevant behaviors and policy support among high- and lower-income individuals.

library(rlang)
library(knitr)
library(Hmisc)
library(dplyr)
library(ggplot2)
library(corrplot)
library(tidyverse)
library(doParallel)
library(BayesFactor)
library(RColorBrewer)
source('helpers.R')

df <- readRDS('Data/BPdata.RDS') %>%
  mutate(
    country = factor(country, c('Denmark', 'United States', 'Nigeria', 'India'))
  )
df$climate_concern <- rowMeans(df[,c('cc_worry', 'cc_importance')], na.rm = TRUE)

actions_order <- c(
  # Curtailment
  'Use less energy for heating and cooling your home',
  'Drive fewer kilometers/miles in your car',
  'Take fewer national and international flights',
  'Eat less red meat',
  'Eat less white meat',
  
  # Investment
  'Improve the insulation of your home',
  'Replace older appliances with newer energy efficient models',
  'Install solar panels at home',
  'Purchase an electric vehicle',
  'Move private investments to climate-friendly financial products'
)

Main Figures

Figure 1: Behavioral plasticity across countries

We look at Behavioral plasticity across countries.

varnames <- c(paste0('curtailment_', seq(5)), paste0('investment_', seq(5)))
map_names <- label(df[varnames])

d <- df %>%
  select(country, varnames) %>%
  pivot_longer(-country) %>% 
  na.omit()

d <- d %>% 
  group_by(country, name, value) %>% 
  summarize(count = n()) %>% 
  left_join(
    d %>% group_by(country, name) %>% summarize(total_count = n()),
    by = c('country', 'name')
  ) %>% 
  mutate(
    prop = count / total_count
  ) %>% 
  mutate(
    value_name = factor(
      ifelse(value == 1, '1: Very unlikely', ifelse(value == 5, '5: Very likely', value)),
      levels = rev(c('1: Very unlikely', '2', '3', '4', '5: Very likely'))
    ),
    full_name = factor(
      wrap_label_vectorized(map_names[name]),
      levels = wrap_label_vectorized(actions_order)
    )
  )

cols <- c('#ff7f0e', '#A65628', '#029E73', '#56B4E9', '#0072B2')
p_country <- ggplot(d, aes(x = country, y = prop, fill = as.factor(value_name))) +
  geom_bar(stat = 'identity') +
  facet_wrap(~ full_name, ncol = 5) + 
  xlab('') +
  ylab('Proportion') +
  scale_fill_manual(
    values = rev(cols)
    # labels = rev(c('1: Very unlikely', '2', '3', '4', '5: Very likely'))
  ) +
  theme_minimal() +
  theme(
    legend.position = 'top',
    panel.grid = element_blank(),
    legend.title = element_blank(),
    strip.text = element_text(size = 7),
    axis.text.x = element_text(angle = 90, hjust = 1),
    plot.title = element_text(size = 14, hjust = 0.50)
  ) +
  guides(fill = guide_legend(reverse = TRUE))

ggsave('Figures/Figure1.pdf', p_country, width = 9, height = 7)
p_country

Calculate proportions for describing the figure.

library(pander)

d_comp <- d %>% 
  group_by(country, full_name) %>% 
  summarize(
    unlikely = round(sum(prop[value %in% c(1, 2)]), 2),
    likely = round(sum(prop[value %in% c(4, 5)]), 2)
  ) %>% data.frame

pander(d_comp)

country	full_name	unlikely	likely
Denmark	Use less energy for heating and cooling your home	0.29	0.45
Denmark	Drive fewer kilometers/miles in your car	0.48	0.27
Denmark	Take fewer national and international flights	0.43	0.33
Denmark	Eat less red meat	0.43	0.38
Denmark	Eat less white meat	0.57	0.2
Denmark	Improve the insulation of your home	0.42	0.36
Denmark	Replace older appliances with newer energy efficient models	0.22	0.53
Denmark	Install solar panels at home	0.57	0.25
Denmark	Purchase an electric vehicle	0.48	0.31
Denmark	Move private investments to climate-friendly financial products	0.46	0.26
United States	Use less energy for heating and cooling your home	0.27	0.48
United States	Drive fewer kilometers/miles in your car	0.36	0.43
United States	Take fewer national and international flights	0.35	0.43
United States	Eat less red meat	0.47	0.33
United States	Eat less white meat	0.59	0.21
United States	Improve the insulation of your home	0.33	0.47
United States	Replace older appliances with newer energy efficient models	0.25	0.49
United States	Install solar panels at home	0.47	0.37
United States	Purchase an electric vehicle	0.5	0.33
United States	Move private investments to climate-friendly financial products	0.45	0.27
Nigeria	Use less energy for heating and cooling your home	0.22	0.58
Nigeria	Drive fewer kilometers/miles in your car	0.35	0.44
Nigeria	Take fewer national and international flights	0.3	0.46
Nigeria	Eat less red meat	0.28	0.53
Nigeria	Eat less white meat	0.44	0.38
Nigeria	Improve the insulation of your home	0.14	0.63
Nigeria	Replace older appliances with newer energy efficient models	0.09	0.73
Nigeria	Install solar panels at home	0.08	0.82
Nigeria	Purchase an electric vehicle	0.3	0.52
Nigeria	Move private investments to climate-friendly financial products	0.16	0.63
India	Use less energy for heating and cooling your home	0.07	0.78
India	Drive fewer kilometers/miles in your car	0.14	0.7
India	Take fewer national and international flights	0.16	0.64
India	Eat less red meat	0.2	0.65
India	Eat less white meat	0.21	0.61
India	Improve the insulation of your home	0.09	0.74
India	Replace older appliances with newer energy efficient models	0.08	0.77
India	Install solar panels at home	0.08	0.82
India	Purchase an electric vehicle	0.09	0.79
India	Move private investments to climate-friendly financial products	0.09	0.73

We use the BayesFactor R package for statistical inference.

Figure 2: Behavioral plasticity and income

We show the relationship between behavioral plasticity and income across items and countries.

df_preds <- do.call('rbind', lapply(varnames, function(varname) get_df_pred_bf(varname, df))) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(map_names[name]),
      levels = wrap_label_vectorized(actions_order)
    )
  )

country_colors <- c('#DE425B', '#3C3B6E', '#008000', '#FF9933')
p_bp_income <- ggplot(df_preds, aes(y = ypred, x = combined_income, fill = country)) +
  geom_vline(aes(xintercept = 10), linetype = 'dotted', color = 'gray48') +
  geom_point(aes(y = mean_value, x = combined_income, col = country), size = 1, show.legend = FALSE) +
  geom_line(aes(col = country), linewidth = 1) +
  geom_ribbon(aes(ymin = ypred_lo, ymax = ypred_hi), alpha = 0.40) + 
  xlab('Income group') +
  ylab('Perceived behavioral plasticity') +
  scale_x_continuous(breaks = c(seq(1, 15, 3), 15), limits = c(1, 15)) +
  scale_y_continuous(breaks = seq(1, 5), limits = c(1, 5)) +
  scale_color_manual(values = country_colors) +
  scale_fill_manual(values = country_colors) +
  facet_wrap(~ full_name) +
  theme_minimal() +
  facet_wrap(~ full_name, ncol = 5) +
  theme(
    plot.title = element_text(size = 14, hjust = 0.50),
    strip.text = element_text(size = 8),
    legend.title = element_blank(),
    legend.position = 'top'
  )

ggsave('Figures/Figure2.pdf', p_bp_income, width = 10, height = 6)
p_bp_income

We calculate inclusion Bayes factors in favor of an effect.

df_inclusion_bfs <- get_bf_inclusion_fig2_df(actions_order)
pander(df_inclusion_bfs %>% select(-rscaleFixed, -rscaleCont))

income	income_x_country	full_name
0.237	0.243	Use less energy for heating and cooling your home
0.396	3.601	Drive fewer kilometers/miles in your car
0.42	5654	Take fewer national and international flights
9.363e+10	0.276	Eat less red meat
2226095	865	Eat less white meat
0.267	4.643	Improve the insulation of your home
7.62	17.21	Replace older appliances with newer energy efficient models
0.154	0.096	Install solar panels at home
3.938e+17	39.3	Purchase an electric vehicle
422.4	0.108	Move private investments to climate-friendly financial products

Figure 3: Behavioral plasticity and domain-matched policy support

We run the relevant three-way interaction models.

df_meat_narm <- df %>%
  filter(!is.na(support_7), !is.na(curtailment_2)) %>% 
  mutate(
    curtailment_2 = curtailment_2 - mean(curtailment_2),
    combined_income = combined_income - mean(combined_income)
  )

df_flights_narm <- df %>%
  filter(!is.na(support_8), !is.na(curtailment_1)) %>% 
  mutate(
    curtailment_1 = curtailment_1 - mean(curtailment_1),
    combined_income = combined_income - mean(combined_income)
  )
  
df_finance_narm <- df %>%
  filter(!is.na(support_6), !is.na(investment_4)) %>% 
  mutate(
    investment_4 = investment_4 - mean(investment_4),
    combined_income = combined_income - mean(combined_income)
  )

rscaleFixed <- sqrt(1/2)
rscaleCont <- sqrt(1/2) / 4

# Lose 422 rows
fit_meat <- generalTestBF(
  support_7 ~ curtailment_2 * combined_income * country,
  data = df_meat_narm, whichModels = 'withmain', 
  rscaleFixed = rscaleFixed, rscaleCont = rscaleCont, progress = FALSE
)

# Lose 443 rows
fit_flights <- generalTestBF(
  support_8 ~ curtailment_1 * combined_income * country,
  data = df_flights_narm, whichModels = 'withmain',
  rscaleFixed = rscaleFixed, rscaleCont = rscaleCont, progress = FALSE
)

# Lose 530 rows
fit_finance <- generalTestBF(
  support_6 ~ investment_4 * combined_income * country,
  data = df_finance_narm, whichModels = 'withmain',
  rscaleFixed = rscaleFixed, rscaleCont = rscaleCont, progress = FALSE
)

# Get the model that predicts the data best
fit_meat_best <- fit_meat[which.max(fit_meat@bayesFactor[, 1])]
fit_flights_best <- fit_flights[which.max(fit_flights@bayesFactor[, 1])]
fit_finance_best <- fit_finance[which.max(fit_finance@bayesFactor[, 1])]

## Used below
countries <- c('Denmark', 'United States', 'Nigeria', 'India')
newdat <- expand.grid(
  country = countries,
  combined_income = seq(15) - mean(seq(15)),
  bp_value = seq(5) - mean(seq(5))
)

map_back <- function(x) {
  as.numeric(factor(x, labels = seq(5)))
}

We look at model predictions using the models with the highest marginal likelihood.

matches <- list(
  c('support_8', 'curtailment_1'),
  c('support_7', 'curtailment_2'),
  c('support_6', 'investment_4')
)

df_preds <- bind_rows(
  get_df_pred_policy(fit_meat_best, 'curtailment_2', 'support_7', df_meat_narm),
  get_df_pred_policy(fit_flights_best, 'curtailment_1', 'support_8', df_flights_narm),
  get_df_pred_policy(fit_finance_best, 'investment_4', 'support_6', df_finance_narm)
) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(map_names[name]),
      levels = wrap_label_vectorized(map_names)
    ),
    
    country = factor(country, levels = c('Denmark', 'United States', 'Nigeria', 'India'))
  )

matching_colors <- c("#E69F00", "#009E73", "#56B4E9")

p_matching <- ggplot(df_preds, aes(x = bp_value, y = ypred, fill = income_bracket)) +
  geom_jitter(
    aes(x = bp_value, y = mean_value, col = income_bracket),
    size = 1, alpha = 1, width = 0.10, show.legend = FALSE
  ) +
  geom_line(aes(col = income_bracket), linewidth = 1) +
  geom_ribbon(aes(ymin = ypred_lo, ymax = ypred_hi), alpha = 0.40, show.legend = FALSE) + 
  xlab('Perceived behavioral plasticity') +
  ylab('Policy support') +
  facet_wrap(~ full_name + country, ncol = 4) +
  scale_y_continuous(breaks = seq(1, 7, 1)) +
  scale_color_manual(values = matching_colors) +
  scale_fill_manual(values = matching_colors) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, hjust = 0.50),
    strip.text = element_text(margin = margin(0.1, 0.1, 0.1, 0.1, 'cm')),
    panel.spacing = unit(0.5, 'lines'),
    legend.position = 'top'
  ) +
  guides(color = guide_legend(title = 'Income bracket'))

ggsave('Figures/Figure3.pdf', p_matching, width = 9, height = 8)
p_matching

library(VGAM)

model_size_prior <- dbetabinom.ab(seq(0, 7), size = 7, shape1 = 1, shape2 = 1)
x <- model_size_prior[1]

# We specify the priors here according to the sequence from the output of the BayesFactor package
# This is if we wanted to use the betabinomial prior, but for the paper we report a uniform model prior
# There's barely any difference
model_prior <- c(
  # Size 0: 1 --> grand mean only
  # Size 1: 3 --> one main effect
  # Size 2: 3 --> two main effects
  # Size 3: 4 --> all main effects or two main effects and one interaction
  # Size 4: 3 --> all main effects and one interactions
  # Size 5: 3 --> all main effects and two two-way interactions
  # Size 6: 1 --> all main effects and all two-way interactions
  # Size 7: 1 --> all main effects and all interactions
  x,
  rep(x / 3, 3),
  rep(x / 3, 3),
  rep(x / 4, 4),
  rep(x / 3, 3),
  rep(x / 3, 3),
  x,
  x
)

set.seed(1)
df_inclusion_bfs3 <- data.frame(rbind(
  get_bf_inclusion_fig3(fit_flights, 'curtailment_1'),
  get_bf_inclusion_fig3(fit_meat, 'curtailment_2'),
  get_bf_inclusion_fig3(fit_finance, 'investment_4')
))

colnames(df_inclusion_bfs3) <- c('bp', 'bp_country', 'bp_income', 'bp_income_country')
df_inclusion_bfs3$full_name <- map_names[c('curtailment_1', 'curtailment_2', 'investment_4')]

pander(df_inclusion_bfs3)

Table continues below

bp	bp_country	bp_income	bp_income_country
2.011e+88	651.9	10.94	0.003961
9.906e+149	123.9	4.204	0.01475
1.36e+101	6800	6124	0.004789

full_name
Take fewer national and international flights
Eat less red meat
Move private investments to climate-friendly financial products

Not so different with uniform prior, which we report in the main paper.

set.seed(1)
df_inclusion_bfs3 <- data.frame(rbind(
  get_bf_inclusion_fig3(fit_flights, 'curtailment_1', prior = 'uniform'),
  get_bf_inclusion_fig3(fit_meat, 'curtailment_2', prior = 'uniform'),
  get_bf_inclusion_fig3(fit_finance, 'investment_4', prior = 'uniform')
))
colnames(df_inclusion_bfs3) <- c('bp', 'bp_country', 'bp_income', 'bp_income_country')
df_inclusion_bfs3$full_name <- map_names[c('curtailment_1', 'curtailment_2', 'investment_4')]

pander(df_inclusion_bfs3)

Table continues below

bp	bp_country	bp_income	bp_income_country
1.661e+88	681.3	11.43	0.003961
8.183e+149	129.5	4.394	0.01475
1.124e+101	7107	6400	0.004789

full_name
Take fewer national and international flights
Eat less red meat
Move private investments to climate-friendly financial products

Extended Data Figures

Figure ED1: Already did the action or cannot do it

We visualize the proportion of people who already did a particular action or cannot do it.

varnames_peb <- paste0(varnames, '_peb')
varnames_cannot <- paste0(varnames, '_cannot')
map_names_cannot <- map_names
names(map_names_cannot) <- varnames_cannot

df_cannot <- df %>%
  select(country, all_of(varnames_cannot)) %>%
  pivot_longer(-country) %>% 
  na.omit() %>% 
  group_by(country, name) %>% 
  summarize(
    prop = mean(value == 1), # Only 0 / 1 exists
    prop_se = sqrt(prop * (1 - prop) / n())
  ) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(map_names_cannot[name]),
      levels = wrap_label_vectorized(actions_order)
    ),
    prop = ifelse(grepl('curtailment', name), NA, prop),
    prop_se = ifelse(grepl('curtailment', name), NA, prop_se)
  )

varnames_peb <- paste0(varnames, '_peb')
map_names_peb <- map_names
names(map_names_peb) <- varnames_peb

df_peb <- df %>%
  select(country, all_of(varnames_peb)) %>%
  pivot_longer(-country) %>% 
  filter(str_starts(name, 'curtailment')) %>% 
  mutate(value = as.numeric(as.character(value))) %>% 
  group_by(country, name) %>% 
  summarize(
    prop = mean(value == -1), 
    prop_se = sqrt(prop * (1 - prop) / n())
  ) %>% 
  bind_rows(
    df %>%
      select(country, all_of(varnames_peb)) %>%
      pivot_longer(-country) %>% 
      mutate(value = as.numeric(as.character(value))) %>% 
      filter(str_starts(name, 'investment')) %>% 
      group_by(country, name) %>% 
      summarize(
        prop = mean(value == 1), 
        prop_se = sqrt(prop * (1 - prop) / n())
      )
  ) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(map_names_peb[name]),
      levels = wrap_label_vectorized(actions_order)
    )
  )

df_both <- bind_rows(
  df_cannot %>% mutate(type = 'Cannot do'),
  df_peb %>% mutate(type = 'Already did / Never do it')
) %>% 
  mutate(
    
    # We change the labels
    full_name = case_when(
      full_name == 'Use less energy for\nheating and cooling your\nhome' ~ 'Do not heat or cool their home',
      full_name == 'Drive fewer kilometers/miles in\nyour car' ~ 'Do not own a car',
      full_name == 'Eat less red meat' ~ 'Do not eat red meat',
      full_name == 'Eat less white meat' ~ 'Do not eat white meat',
      full_name == 'Take fewer national and\ninternational flights' ~ 'Never fly',
      TRUE ~ full_name
    ),
    
    full_name = factor(
      full_name,
      levels = c(
        rev(unname(wrap_label_vectorized(actions_order)))[seq(5)],
        'Do not eat white meat',
        'Do not eat red meat',
        'Never fly',
        'Do not own a car',
        'Do not heat or cool their home'
      )
    )
  )

pwidth <- 0.96
p_both <- ggplot(df_both, aes(x = prop, y = full_name, fill = country, col = country)) +
  geom_bar(stat = 'identity', position = position_dodge(width = pwidth)) +
    geom_point(
    aes(x = prop, y = full_name), position = position_dodge(width = pwidth),
    size = 2, show.legend = FALSE, color = 'black'
  ) +
  geom_errorbar(
    aes(xmin = prop - 1.96 * prop_se, xmax = prop + 1.96 * prop_se),
    position = position_dodge(width = pwidth),
    width = 0.40, linewidth = 1, show.legend = FALSE, color = 'black'
  ) +
  geom_point(
    aes(x = prop, y = full_name), position = position_dodge(width = pwidth),
    size = 1, show.legend = FALSE
  ) +
  geom_errorbar(
    aes(xmin = prop - 1.96 * prop_se, xmax = prop + 1.96 * prop_se),
    position = position_dodge(width = pwidth),
    width = 0.30, linewidth = 0.30,
    show.legend = FALSE
  ) +
  facet_wrap(~ type) +
  ylab('') +
  xlab('Proportion') +
  scale_fill_manual(values = country_colors) +
  scale_color_manual(values = country_colors) +
  theme_minimal() +
  theme(
    legend.position = 'top',
    legend.title = element_blank(),
    strip.text = element_text(size = 10),
    plot.title = element_text(size = 14, hjust = 0.50)
  )

ggsave('Figures/FigureED1.pdf', p_both, width = 9, height = 8)
p_both

Figure ED2: Already did and behavioral plasticity

Visualize relationship between exclusion rate and mean plasticity rating for behaviors across countries.

df_mean <- df %>%
  select(country, varnames) %>%
  pivot_longer(-country) %>% 
  na.omit() %>% 
  group_by(country, name) %>% 
  summarize(
    mean_bp = mean(value),
    sd_bp = sd(value),
    se_bp = sd_bp / sqrt(n())
  ) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(map_names[name]),
      levels = wrap_label_vectorized(actions_order)
    )
  )

df_mean_exclusion <- df_mean %>% 
  left_join(
    df_both %>%
      filter(type == 'Already did / Never do it') %>% 
      mutate(
        name = str_extract(name, "^[^_]+_[^_]+")
      ) %>% 
      select(country, name, prop, prop_se)
    , by = c('country', 'name')
  ) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(map_names[name]),
      levels = rev(wrap_label_vectorized(actions_order))
    )
  )
  
cols <- brewer.pal(10, 'Set3')
p_mean_exclusion <- ggplot(df_mean_exclusion, aes(x = prop, y = mean_bp, col = full_name)) +
  geom_point(size = 2) +
  geom_errorbar(
    aes(ymin = mean_bp - 1.96 * se_bp, ymax = mean_bp + 1.96 * se_bp),
    width = 0.01
  ) +
  # Horizontal error bars
  geom_errorbarh(
    aes(xmin = prop - 1.96 * prop_se, xmax = prop + 1.96 * prop_se),
    height = 0.1
  ) +
  facet_wrap(~ country) +
  labs(
    x = 'Proportion: Already do / never do it',
    y = 'Mean perceived behavioral plasticity'
  ) +
  scale_color_manual(name = '', values = cols) +
  scale_x_continuous(breaks = seq(0, 0.3, 0.05)) +
  scale_y_continuous(breaks = seq(1, 5, 1), limits = c(1, 5)) +
  theme_minimal() +
  theme(
    legend.position = 'top',
    strip.text = element_text(size = 14),
    plot.title = element_text(hjust = 0.50, size = 14),
    axis.text.y = element_text(size = 10),
    legend.key.spacing = unit(0.10, 'cm')
  ) +
  guides(color = guide_legend(nrow = 5, reverse = TRUE))

ggsave('Figures/FigureED2.pdf', p_mean_exclusion, width = 7, height = 7)
p_mean_exclusion

Figure ED3: Already did the action across income groups

df_prep <- df %>%
  select(country, combined_income, varnames_peb) %>%
  pivot_longer(cols = -c(country, combined_income)) %>% 
  mutate(value = as.numeric(as.character(value))) %>% 
  mutate(
    income_bracket = ifelse(
      combined_income <= 5, 'low',
      ifelse(combined_income < 10, 'medium', 'high')
    )
  )

df_selection <- df_prep %>% 
  filter(str_starts(name, 'curtailment')) %>% 
  group_by(country, income_bracket, name) %>% 
  summarize(
    prop = mean(value == -1), 
    prop_se = sqrt(prop * (1 - prop) / n())
  ) %>% 
  bind_rows(
    df_prep %>%
      filter(str_starts(name, 'investment')) %>% 
      group_by(country, income_bracket, name) %>% 
      summarize(
        prop = mean(value == 1), 
        prop_se = sqrt(prop * (1 - prop) / n())
      )
  ) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(map_names_peb[name]),
      levels = wrap_label_vectorized(actions_order)
    ),
    
    income_bracket = factor(
      income_bracket,
      levels = c('low', 'medium', 'high'),
      labels = c('Low (1 - 5)', 'Medium (6 - 9)', 'High (10 - 15)')
    ),
    
    # We change the labels
    full_name = case_when(
      full_name == 'Use less energy for\nheating and cooling your\nhome' ~ 'Do not heat or cool their home',
      full_name == 'Drive fewer kilometers/miles in\nyour car' ~ 'Do not own a car',
      full_name == 'Eat less red meat' ~ 'Do not eat red meat',
      full_name == 'Eat less white meat' ~ 'Do not eat white meat',
      full_name == 'Take fewer national and\ninternational flights' ~ 'Never fly',
      TRUE ~ full_name
    ),
    
    full_name = factor(
      full_name,
      levels = c(
        'Do not heat or cool their home',
        'Do not own a car',
        'Never fly',
        'Do not eat red meat',
        'Do not eat white meat',
        rev(unname(wrap_label_vectorized(actions_order)))[seq(5)]
      )
    )
  )

p_selection <- ggplot(
    df_selection,
    aes(y = prop, x = country, fill = income_bracket)
  ) +
  geom_bar(stat = 'identity', position = position_dodge(width = pwidth)) +
  geom_point(
    aes(x = country, y = prop), position = position_dodge(width = pwidth),
    size = 2, show.legend = FALSE, color = 'black'
  ) +
  geom_errorbar(
    aes(ymin = prop - 1.96 * prop_se, ymax = prop + 1.96 * prop_se),
    position = position_dodge(width = pwidth),
    width = 0.40, linewidth = 1, show.legend = FALSE, color = 'black'
  ) +
  geom_point(
    aes(x = country, y = prop, col = income_bracket), position = position_dodge(width = pwidth),
    size = 1, show.legend = FALSE
  ) +
  geom_errorbar(
    aes(ymin = prop - 1.96 * prop_se, ymax = prop + 1.96 * prop_se, col = income_bracket),
    position = position_dodge(width = pwidth),
    width = 0.30, linewidth = 0.30,
    show.legend = FALSE
  ) +
  facet_wrap(~ full_name, nrow = 2) +
  ylab('Proportion') +
  xlab('') +
  scale_fill_manual(values = matching_colors) +
  scale_color_manual(values = matching_colors) +
  theme_minimal() +
  theme(
    legend.position = 'top',
    legend.title = element_blank(),
    strip.text = element_text(size = 7),
    axis.text.x = element_text(angle = 90),
    plot.title = element_text(size = 14, hjust = 0.50)
  )

ggsave('Figures/FigureEx3.pdf', p_selection, width = 9, height = 7)
p_selection

Let’s run a proportion test to assess the evidence for differences.

df_bfs <- c()
for (country in unique(df_prep$country)) {
  for (name in unique(df_prep$name)) {
    bf <- test_prop(df, country, name)
    df_bfs <- rbind(df_bfs, cbind(bf, country, name))
  }
}

df_bfs <- data.frame(df_bfs) %>% 
  mutate(
    name = gsub('_peb$', '', name),
    log_bf = as.numeric(bf),
    full_name = map_names[name]
  )

pander(df_bfs %>% select(country, full_name, log_bf))

country	full_name	log_bf
Denmark	Take fewer national and international flights	45.66
Denmark	Eat less red meat	-5.937
Denmark	Eat less white meat	-6.544
Denmark	Drive fewer kilometers/miles in your car	24.79
Denmark	Use less energy for heating and cooling your home	-5.283
Denmark	Purchase an electric vehicle	25.19
Denmark	Replace older appliances with newer energy efficient models	3.275
Denmark	Improve the insulation of your home	22.6
Denmark	Move private investments to climate-friendly financial products	-1.516
Denmark	Install solar panels at home	-0.9868
United States	Take fewer national and international flights	52.32
United States	Eat less red meat	-5.661
United States	Eat less white meat	-6.602
United States	Drive fewer kilometers/miles in your car	17.95
United States	Use less energy for heating and cooling your home	-6.944
United States	Purchase an electric vehicle	-0.01785
United States	Replace older appliances with newer energy efficient models	0.4505
United States	Improve the insulation of your home	6.468
United States	Move private investments to climate-friendly financial products	-5.966
United States	Install solar panels at home	0.8944
India	Take fewer national and international flights	-1.688
India	Eat less red meat	-4.436
India	Eat less white meat	-4.312
India	Drive fewer kilometers/miles in your car	-1.241
India	Use less energy for heating and cooling your home	-5.221
India	Purchase an electric vehicle	1.757
India	Replace older appliances with newer energy efficient models	-4.226
India	Improve the insulation of your home	-3.277
India	Move private investments to climate-friendly financial products	-4.846
India	Install solar panels at home	-4.477
Nigeria	Take fewer national and international flights	7.376
Nigeria	Eat less red meat	-5.844
Nigeria	Eat less white meat	3.438
Nigeria	Drive fewer kilometers/miles in your car	8.633
Nigeria	Use less energy for heating and cooling your home	-1.016
Nigeria	Purchase an electric vehicle	-4.526
Nigeria	Replace older appliances with newer energy efficient models	25.45
Nigeria	Improve the insulation of your home	-4.727
Nigeria	Move private investments to climate-friendly financial products	-5.887
Nigeria	Install solar panels at home	30.04

Supplementary Figures

Figure S1: Mean behavioral plasticity across countries

Here we visualize mean behavioral plasticity across all countries.

df_mean_sd_beh <- df %>%
  select(country, varnames) %>%
  pivot_longer(-country) %>% 
  na.omit() %>% 
  mutate(curtailment = as.numeric(grepl('curtailment', name))) %>% 
  group_by(country, name) %>% 
  summarize(
    mean_bp = mean(value),
    sd_bp = sd(value),
    n = n(),
    se_bp = sd_bp / sqrt(n)
  ) %>% 
  # Compute the mean of the mean /sd behavioral plasticity of behaviors
  group_by(country) %>% 
  summarize(
    mean_bp = mean(mean_bp),
    sd_bp = mean(sd_bp)
  )

pander(df_mean_sd_beh)

country	mean_bp	sd_bp
Denmark	2.777	1.365
United States	2.902	1.461
Nigeria	3.563	1.292
India	4.024	1.149

df_mean <- df %>%
  select(country, varnames) %>%
  pivot_longer(-country) %>% 
  na.omit() %>% 
  group_by(country, name) %>% 
  summarize(
    mean_bp = mean(value),
    sd_bp = sd(value),
    n = n(),
    se_bp = sd_bp / sqrt(n),
    se_sd_bp = sd_bp / sqrt(2 * n)
  ) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(map_names[name]),
      levels = wrap_label_vectorized(actions_order)
    )
  )

pwidth <- 0.91
p_mean <- ggplot(df_mean, aes(x = mean_bp, y = country, color = country)) +
  geom_point(position = position_dodge(width = pwidth)) +
  geom_errorbar(
    aes(xmin = mean_bp - 1.96 * se_bp, xmax = mean_bp + 1.96 * se_bp),
    position = position_dodge(width = pwidth), width = 0.30
  ) +
  facet_wrap(~ full_name, nrow = 2) +
  ylab('') +
  xlab('Mean perceived behavioral plasticity') +
  scale_x_continuous(breaks = seq(1, 5, 1), limits = c(1, 5)) +
  scale_color_manual(values = country_colors) +
  coord_flip() +
  theme_minimal() +
  theme(
    legend.position = 'top',
    legend.title = element_blank(),
    axis.text.x = element_text(angle = 90),
    strip.text = element_text(size = 7),
    panel.grid.major.x = element_blank(),
    plot.title = element_text(size = 14, hjust = 0.50)
  )

ggsave('Figures/FigureS1.pdf', p_mean, width = 9, height = 7)
p_mean

Figure S1: Standard deviation of behavioral plasticity across countries

pwidth <- 0.91
p_sd <- ggplot(df_mean, aes(x = sd_bp, y = country, color = country)) +
  geom_point(position = position_dodge(width = pwidth)) +
  geom_errorbar(
    aes(xmin = sd_bp - 1.96 * se_sd_bp, xmax = sd_bp + 1.96 * se_sd_bp),
    position = position_dodge(width = pwidth), width = 0.30
  ) +
  facet_wrap(~ full_name, nrow = 2) +
  ylab('') +
  xlab('Standard deviation of perceived behavioral plasticity') +
  scale_x_continuous(breaks = seq(0, 2, 0.50), limits = c(0, 2)) +
  scale_color_manual(values = country_colors) +
  coord_flip() +
  theme_minimal() +
  theme(
    legend.position = 'top',
    legend.title = element_blank(),
    axis.text.x = element_text(angle = 90),
    strip.text = element_text(size = 7),
    panel.grid.major.x = element_blank(),
    plot.title = element_text(size = 14, hjust = 0.50)
  )

ggsave('Figures/FigureS2.pdf', p_sd, width = 9, height = 7)
p_sd

Figure S3: Mean policy support across countries

Here we visualize mean policy support across all countries.

support_vars <- colnames(df %>% select(starts_with('support_')))
support_labels <- label(df %>% select(starts_with('support_')))

df_mean_ps <- df %>%
  select(country, starts_with('support_')) %>%
  pivot_longer(-country) %>% 
  na.omit() %>% 
  group_by(country, name) %>% 
  summarize(
    mean_ps = mean(value),
    sd_ps = sd(value),
    n = n(),
    se_ps = sd_ps / sqrt(n)
  ) %>% 
  mutate(
    full_name = factor(
      wrap_label_vectorized(support_labels[name]),
      levels = wrap_label_vectorized(support_labels[name])
    )
  )

pwidth <- 0.91
country_colors <- c('#DE425B', '#3C3B6E', '#008000', '#FF9933')
p_mean_ps <- ggplot(df_mean_ps, aes(x = mean_ps, y = country, color = country)) +
  geom_point(position = position_dodge(width = pwidth)) +
  geom_errorbar(
    aes(xmin = mean_ps - 1.96 * se_ps, xmax = mean_ps + 1.96 * se_ps),
    position = position_dodge(width = pwidth), width = 0.30
  ) +
  facet_wrap(~ full_name, nrow = 2) +
  ylab('') +
  xlab('Mean policy support') +
  scale_x_continuous(breaks = seq(1, 7, 1), limits = c(1, 7)) +
  scale_color_manual(values = country_colors) +
  coord_flip() +
  theme_minimal() +
  theme(
    legend.position = 'top',
    legend.title = element_blank(),
    axis.text.x = element_text(angle = 90),
    strip.text = element_text(size = 6),
    panel.grid.major.x = element_blank(),
    plot.title = element_text(size = 14, hjust = 0.50)
  )

ggsave('Figures/FigureS3.pdf', p_mean_ps, width = 9, height = 8)
p_mean_ps

Figure S4: Differences in high and low income groups

Here we visualize the differences in mean behavioral plasticity between the top 10% income earners and the rest.

set.seed(1)
df_diffs <- do.call(
  'rbind', lapply(varnames, function(varname) get_top10_difference(varname, df))
) %>% 
  mutate(
    country = factor(country, levels = rev(countries)),
    full_name = factor(
      map_names[name],
      levels = rev(actions_order)
    )
  )

cols <- brewer.pal(10, 'Set3')
p_estimates <- ggplot(df_diffs, aes(x = mean_diff, y = country, color = full_name)) +
  geom_vline(aes(xintercept = 0), linetype = 'dotted', color = 'gray48') +
  geom_errorbar(
    aes(xmin = mean_diff_lo, xmax = mean_diff_hi),
    position = position_dodge(width = 0.70), width = 0.80, size = 1.5
  ) +
  geom_point(position = position_dodge(width = 0.70), size = 3) + 
  scale_x_continuous(breaks = round(seq(-0.80, 1.2, 0.2), 3), limits = c(-0.90, 1.3)) +
  theme_minimal() +
  scale_color_manual(name = '', values = cols) +
  labs(
    x = 'Mean difference in perceived behavioral plasticity',
    y = ''
  ) +
  theme(
    legend.position = 'top',
    strip.text = element_text(size = 14),
    plot.title = element_text(hjust = 0.50, size = 14),
    axis.text.y = element_text(size = 10),
    legend.key.spacing = unit(0.06, 'cm')
  ) +
  guides(color = guide_legend(nrow = 5, reverse = TRUE))

ggsave('Figures/FigureS4.pdf', p_estimates, width = 8, height = 10)
p_estimates

Figures S5 & S6: Sensitivity analysis I

We assess how sensitive the results are to changes in the prior width of the income effect and the interaction of income and country.

scale_lo <- 0.05
scale_hi <- 0.30
scales_cont <- seq(scale_lo, scale_hi, 0.0125)

df_inclusion_bfs_sens <- do.call('rbind', lapply(scales_cont, function(scale) {
  get_bf_inclusion_fig2_df(actions_order, rscaleCont = scale)
}))
end <- Sys.time()

df_inclusion_bfs_sens <- df_inclusion_bfs_sens %>% 
    mutate(
    full_name = factor(
      wrap_label_vectorized(full_name),
      levels = wrap_label_vectorized(actions_order)
    )
  )

library(scales)

# Custom function to apply scientific notation for numbers >= 10
custom_label <- function(x) {
  ifelse(x > 100, scientific_format()(x), format(x, scientific = FALSE))
}

df_inclusion_bfs_sens2 <- df_inclusion_bfs_sens %>% 
  mutate(income = ifelse(income <= 1, 1 / income, income))

p_sens_fig2_income <- ggplot(df_inclusion_bfs_sens, aes(x = rscaleCont, y = (income))) +
  geom_vline(aes(xintercept = sqrt(1/2) / 4), linetype = 'dotted', color = 'gray48') +
  labs(
    y = expression('Inclusion ' ~ BF[1*0]),
    x = 'Cauchy prior width'
  ) +
  geom_line(linewidth = 1) +
  scale_y_continuous(labels = custom_label) +
  scale_x_continuous(breaks = seq(0.05, 0.30, 0.05)) +
  facet_wrap(~ full_name, ncol = 5, scales = 'free_y') +
  theme_minimal() +
  theme(
    legend.position = 'none',
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    strip.text = element_text(size = 7)
  )

ggsave(
  'Figures/FigureS5.pdf',
  p_sens_fig2_income,
  width = 10, height = 5
)
p_sens_fig2_income

p_sens_fig2_income_country <- ggplot(
  df_inclusion_bfs_sens, aes(x = rscaleCont, y = (income_x_country))
) +
  geom_vline(aes(xintercept = sqrt(1/2) / 4), linetype = 'dotted', color = 'gray48') +
  labs(
    y = expression('Inclusion ' ~ BF[1*0]),
    x = 'Cauchy prior width'
  ) +
  geom_line(linewidth = 1) +
  scale_y_continuous(labels = custom_label) +
  scale_x_continuous(breaks = seq(0.05, 0.30, 0.05)) +
  facet_wrap(~ full_name, ncol = 5, scales = 'free_y') +
  theme_minimal() +
  theme(
    legend.position = 'none',
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    strip.text = element_text(size = 7)
  )

ggsave(
  'Figures/FigureS6.pdf',
  p_sens_fig2_income_country,
  width = 10, height = 5
)
p_sens_fig2_income_country

Figure S7: Sensitivity analysis II

We run sensitivity analyses for the results reported in Figure 3.

if (!file.exists('Data/df_sensitivity_matched.csv')) {
  
  rscaleFixed <- sqrt(1/2)
  scales_cont <- seq(0.05, 0.30, 0.0125)
  
  df_inclusion_bfs <- do.call('rbind', lapply(scales_cont, function(scale) {
    
    set.seed(1)
    fit_meat <- generalTestBF(
      support_7 ~ curtailment_2 * combined_income * country,
      data = df_meat_narm, whichModels = 'withmain', 
      rscaleFixed = rscaleFixed, rscaleCont = scale, progress = FALSE
    )
    
    set.seed(1)
    fit_flights <- generalTestBF(
      support_8 ~ curtailment_1 * combined_income * country,
      data = df_flights_narm, whichModels = 'withmain',
      rscaleFixed = rscaleFixed, rscaleCont = scale, progress = FALSE
    )
    
    set.seed(1)
    fit_finance <- generalTestBF(
      support_6 ~ investment_4 * combined_income * country,
      data = df_finance_narm, whichModels = 'withmain',
      rscaleFixed = rscaleFixed, rscaleCont = scale, progress = FALSE
    )
    
    df_inclusion_bfs <- data.frame(rbind(
      get_bf_inclusion_fig3(fit_flights, 'curtailment_1', prior = 'uniform'),
      get_bf_inclusion_fig3(fit_meat, 'curtailment_2', prior = 'uniform'),
      get_bf_inclusion_fig3(fit_finance, 'investment_4', prior = 'uniform')
    ))
    
    colnames(df_inclusion_bfs) <- c('bp', 'bp_country', 'bp_income', 'bp_income_country')
    df_inclusion_bfs$full_name <- map_names[c('curtailment_1', 'curtailment_2', 'investment_4')]
    df_inclusion_bfs$rscaleCont <- scale
    df_inclusion_bfs
    
  }))
  
  df_inclusion_bfs_sens2 <- df_inclusion_bfs %>% 
    mutate(
      full_name = factor(
        wrap_label_vectorized(full_name),
        levels = wrap_label_vectorized(actions_order)
      )
    )
  
  df_inclusion_bfs_sens2_long <- df_inclusion_bfs_sens2 %>%
    pivot_longer(
      cols = c(bp, bp_country, bp_income, bp_income_country),
      names_to = 'type',
      values_to = 'value'
    ) %>% 
    mutate(
      type = case_when(
        type == 'bp' ~ 'BP',
        type == 'bp_country' ~ 'BP x Country',
        type == 'bp_income' ~ 'BP x Income',
        type == 'bp_income_country' ~ 'BP x Income x Country'
      )
    )
} else {
  df_inclusion_bfs_sens2_long <- read.csv('Data/df_sensitivity_matched.csv') %>% 
    mutate(
      full_name = factor(
        full_name, levels = c(
          'Take fewer national and\ninternational flights', 'Eat less red meat',
          'Move private investments to\nclimate-friendly financial products'
        )
      )
    )
}

p_sens_fig3 <- ggplot(
  df_inclusion_bfs_sens2_long, aes(x = rscaleCont, y = value)
) +
  geom_vline(aes(xintercept = sqrt(1/2) / 4), linetype = 'dotted', color = 'gray48') +
  labs(
    y = expression('Inclusion ' ~ BF[1*0]),
    x = 'Cauchy prior width'
  ) +
  geom_line(linewidth = 1) +
  scale_y_continuous(labels = label_scientific()) +
  scale_x_continuous(breaks = seq(0.05, 0.30, 0.05)) +
  # facet_nested(~ type + full_name, scales = 'free_y') +
  facet_wrap(~ type + full_name, ncol = 3, scales = 'free_y') +
  theme_minimal() +
  theme(
    legend.position = 'none',
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    strip.text = element_text(size = 7)
  )

ggsave(
  'Figures/FigureS7.pdf',
  p_sens_fig3,
  width = 8, height = 8
)
p_sens_fig3

Figures S8 & S9: More loosely domain-matched analyses

We estimate the relationship between more loosely domain-matched policies and behavioral plasticity items.

get_best <- function(fit_all) fit_all[which.max(fit_all@bayesFactor[, 1])]

matches_loose <- list(
  
  # Increasing energy efficiency requirements for buildings & insulating your home
  c('support_5', 'investment_3'),
  
  # Increase price of electricity during peak hours & replace older with more energy efficient appliances
  c('support_3', 'investment_2'),
  
  # increasing subsidies for renewable energy projects & install solar panels at home
  c('support_4', 'investment_5'),
  
  # Ban sale of diesel & purchasing an EV
  c('support_10', 'investment_1'),
  
  # Requires banks to disclose & move investments
  c('support_6', 'investment_4'),
  
  # Increase subsidies for food products with low GHGs & eating less red meat
  c('support_12', 'curtailment_2'),
  
  # Increase subsidies for food products with low GHGs & eating less red meat
  c('support_12', 'curtailment_3'),
  
  # Expand public transport & drive fewer kilometres with the car
  c('support_2', 'curtailment_4'),
  
  # Strengthen requirements for energy efficiency & use less energy for heating and cooling
  c('support_5', 'curtailment_5')
)

#' Prepare data.frame for loosely domain-matched policy analysis
prepare_df <- function(df, bp_name, policy_name) {
  bp_name <- ensym(bp_name)
  policy_name <- ensym(policy_name)
  
  df %>% 
    filter(!is.na(!!bp_name), !is.na(!!policy_name)) %>% 
    mutate(
      !!bp_name := !!bp_name - mean(!!bp_name, na.rm = TRUE),
      combined_income = combined_income - mean(combined_income, na.rm = TRUE)
    )
}

if (!file.exists('Data/df_preds_loose.csv')) {
  
  fit_cars <- get_bf_fit_loose(df, 'support_10', 'investment_1')
  fit_electricity <- get_bf_fit_loose(df, 'support_3', 'investment_2')
  fit_banks <- get_bf_fit_loose(df, 'support_6', 'investment_4')
  fit_replace <- get_bf_fit_loose(df, 'support_5', 'investment_3')
  fit_renewables <- get_bf_fit_loose(df, 'support_4', 'investment_5')
  
  fit_food_red <- get_bf_fit_loose(df, 'support_12', 'curtailment_2')
  fit_food_white <- get_bf_fit_loose(df, 'support_12', 'curtailment_3')
  fit_public <- get_bf_fit_loose(df, 'support_2', 'curtailment_4')
  fit_efficiency <- get_bf_fit_loose(df, 'support_5', 'curtailment_5')
  
  df_preds_loose <- bind_rows(
    get_df_pred_policy(
      get_best(fit_cars), 'investment_1', 'support_10',
      prepare_df(df, 'investment_1', 'support_10')
    ) %>% mutate(type = 'cars'),
    
    get_df_pred_policy(
      get_best(fit_electricity), 'investment_2', 'support_3',
      prepare_df(df, 'investment_2', 'support_3')
    ) %>% mutate(type = 'electricity'),
    
    get_df_pred_policy(
      get_best(fit_banks), 'investment_4', 'support_6',
      prepare_df(df, 'investment_4', 'support_6')
    ) %>% mutate(type = 'banks'),
    
    get_df_pred_policy(
      get_best(fit_replace), 'investment_3', 'support_5',
      prepare_df(df, 'investment_3', 'support_5')
    ) %>% mutate(type = 'replace'),
    
    get_df_pred_policy(
      get_best(fit_renewables), 'investment_5', 'support_4',
      prepare_df(df, 'investment_5', 'support_4')
    ) %>% mutate(type = 'renewables'),
    
    get_df_pred_policy(
      get_best(fit_food_red), 'curtailment_2', 'support_12',
      prepare_df(df, 'curtailment_2', 'support_12')
    ) %>% mutate(type = 'red_meat'),
    
    get_df_pred_policy(
      get_best(fit_food_white), 'curtailment_3', 'support_12',
      prepare_df(df, 'curtailment_3', 'support_12')
    ) %>% mutate(type = 'white_meat'),
    
    get_df_pred_policy(
      get_best(fit_public), 'curtailment_4', 'support_2',
      prepare_df(df, 'curtailment_4', 'support_2')
    ) %>% mutate(type = 'public'),
    
    get_df_pred_policy(
      get_best(fit_efficiency), 'curtailment_5', 'support_5',
      prepare_df(df, 'curtailment_5', 'support_5')
    ) %>% mutate(type = 'public')
  ) %>% 
    mutate(
      full_name = factor(
        wrap_label_vectorized(map_names[name]),
        levels = wrap_label_vectorized(actions_order)
      ),
      country = factor(country, levels = c('Denmark', 'United States', 'Nigeria', 'India'))
    )
} else {
  df_preds_loose <- read.csv('Data/df_preds_loose.csv') %>% 
    mutate(
      full_name = factor(
        wrap_label_vectorized(map_names[name]),
        levels = wrap_label_vectorized(actions_order)
      ),
      country = factor(country, levels = c('Denmark', 'United States', 'Nigeria', 'India')),
      income_bracket = factor(income_bracket, levels = c('Low (1 - 5)', 'Medium (6 - 9)', 'High (10 - 15)'))
    )
}

df_preds_curtailment_loose <- df_preds_loose %>% 
  filter(str_starts(name, 'curtailment'))

df_preds_investment_loose <- df_preds_loose %>% 
  filter(str_starts(name, 'investment'))

matching_colors <- c("#E69F00", "#009E73", "#56B4E9")

p_investment_loose <- ggplot(
  df_preds_investment_loose, aes(x = bp_value, y = ypred, fill = income_bracket)) +
  geom_jitter(
    aes(x = bp_value, y = mean_value, col = income_bracket),
    size = 1, alpha = 1, width = 0.10, show.legend = FALSE
  ) +
  geom_line(aes(col = income_bracket), linewidth = 1) +
  geom_ribbon(aes(ymin = ypred_lo, ymax = ypred_hi), alpha = 0.40, show.legend = FALSE) + 
  xlab('Perceived behavioral plasticity') +
  ylab('Policy support') +
  facet_wrap(~ full_name + country, ncol = 4) +
  scale_color_manual(values = matching_colors) +
  scale_fill_manual(values = matching_colors) +
  theme_minimal() +
  scale_y_continuous(breaks = seq(1, 7), limits = c(1, 7)) +
  theme(
    plot.title = element_text(size = 14, hjust = 0.50),
    strip.text = element_text(margin = margin(0.1, 0.1, 0.1, 0.1, 'cm')),
    panel.spacing = unit(0.5, 'lines'),
    legend.position = 'top'
  ) +
  guides(color = guide_legend(title = 'Income bracket'))

p_curtailment_loose <- ggplot(
  df_preds_curtailment_loose, aes(x = bp_value, y = ypred, fill = income_bracket)) +
  geom_jitter(
    aes(x = bp_value, y = mean_value, col = income_bracket),
    size = 1, alpha = 1, width = 0.10, show.legend = FALSE
  ) +
  geom_line(aes(col = income_bracket), linewidth = 1) +
  geom_ribbon(aes(ymin = ypred_lo, ymax = ypred_hi), alpha = 0.40, show.legend = FALSE) + 
  xlab('Behavioral plasticity') +
  ylab('Policy support') +
  scale_y_continuous(breaks = seq(1, 7), limits = c(1, 7)) +
  facet_wrap(~ full_name + country, ncol = 4) +
  scale_color_manual(values = matching_colors) +
  scale_fill_manual(values = matching_colors) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, hjust = 0.50),
    strip.text = element_text(margin = margin(0.1, 0.1, 0.1, 0.1, 'cm')),
    panel.spacing = unit(0.5, 'lines'),
    legend.position = 'top'
  ) +
  guides(color = guide_legend(title = 'Income bracket'))

ggsave('Figures/FigureS8.pdf', p_curtailment_loose, width = 9, height = 10)
p_curtailment_loose

ggsave('Figures/FigureS9.pdf', p_investment_loose, width = 9, height = 10)
p_investment_loose

Figure S10: Mean behavioral plasticy and mean policy support

We run the model for mean support across all policies and mean of all behavioral plasticity items.

df_avg <- df %>% 
  rowwise() %>% 
  mutate(
    support_avg = mean(c_across(starts_with('support')), na.rm = TRUE),
    bp_value = mean(c_across(c(paste0('curtailment_', seq(5)), paste0('investment_', seq(5)))), na.rm = TRUE)
  ) %>% data.frame

fit_all <- generalTestBF(
  support_avg ~ bp_value * combined_income * country,
  data = df_avg, whichModels = 'withmain',
  rscaleFixed = rscaleFixed, rscaleCont = rscaleCont, progress = FALSE
)

fit_all_best <- fit_all[which.max(fit_all@bayesFactor[, 1])]

df_preds_avg <- get_df_pred_policy(
  fit_all_best, 'bp_value', 'support_avg', df_avg, type = 'averaged',
  join = FALSE, use_empirical = FALSE
) %>% 
  mutate(
    country = factor(country, levels = c('Denmark', 'United States', 'Nigeria', 'India')),
    income_bracket = factor(
      income_bracket,
      levels = c('low', 'medium', 'high'),
      labels = c('Low (1 - 5)', 'Medium (6 - 9)', 'High (10 - 15)')
    )
  )

df_emp_avg <- df_avg %>% 
  mutate(
    income_bracket = ifelse(
      combined_income <= 5, 'low',
      ifelse(combined_income < 10, 'medium', 'high')
    ),
    income_bracket = factor(
      income_bracket,
      levels = c('low', 'medium', 'high'),
      labels = c('Low (1 - 5)', 'Medium (6 - 9)', 'High (10 - 15)')
    )
  )

df_both <- df_emp_avg %>% 
  left_join(df_preds_avg, by = c('country', 'bp_value', 'income_bracket'))

p_mean_matched <- ggplot(df_both, aes(x = bp_value, y = ypred, fill = income_bracket)) +
  geom_jitter(
    aes(x = bp_value, y = support_avg, col = income_bracket),
    size = 1, alpha = 0.30, width = 0.10, show.legend = FALSE
  ) +
  geom_line(data = df_preds_avg, aes(col = income_bracket), linewidth = 1) +
  geom_ribbon(data = df_preds_avg, aes(ymin = ypred_lo, ymax = ypred_hi), alpha = 0.40, show.legend = FALSE) + 
  scale_y_continuous(breaks = seq(1, 7, 1), limits = c(0.9, 7.2)) +
  xlab('Mean perceived behavioral plasticity') +
  ylab('Mean policy support') +
  facet_wrap(~ country, ncol = 4) +
  scale_color_manual(values = matching_colors) +
  scale_fill_manual(values = matching_colors) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, hjust = 0.50),
    strip.text = element_text(margin = margin(0.1, 0.1, 0.1, 0.1, 'cm')),
    panel.spacing = unit(0.5, 'lines'),
    legend.position = 'top'
  ) +
  guides(color = guide_legend(title = 'Income bracket'))

ggsave('Figures/FigureS10.pdf', p_mean_matched, width = 9, height = 4)
p_mean_matched

Supplementary Figures 11 - 14: Correlation plots

We calculate Kendall’s correlation coefficients.

support_vars <- colnames(df %>% select(starts_with('support_')))
support_labels <- label(df %>% select(starts_with('support_')))

df_cors <- df %>% select(all_of(varnames), starts_with('support_'), country)
labels <- c(wrap_label_vectorized(map_names), wrap_label_vectorized(support_labels))

res_india <- get_cors_all(df_cors, labels, 'India')
res_denmark <- get_cors_all(df_cors, labels, 'Denmark')
res_nigeria <- get_cors_all(df_cors, labels, 'Nigeria')
res_usa <- get_cors_all(df_cors, labels, 'United States')

Denmark

pdf('Figures/FigureS11.pdf', width = 10, height = 10)
get_corrplot_all(res_denmark)
dev.off()

## quartz_off_screen 
##                 2

get_corrplot_all(res_denmark)

United States

pdf('Figures/FigureS12.pdf', width = 10, height = 10)
get_corrplot_all(res_usa)
dev.off()

## quartz_off_screen 
##                 2

get_corrplot_all(res_usa)

Nigeria

pdf('Figures/FigureS13.pdf', width = 10, height = 10)
get_corrplot_all(res_nigeria)
dev.off()

## quartz_off_screen 
##                 2

get_corrplot_all(res_nigeria)

India

pdf('Figures/FigureS14.pdf', width = 10, height = 10)
get_corrplot_all(res_india)
dev.off()

## quartz_off_screen 
##                 2

get_corrplot_all(res_india)